{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Opencast Mining\n",
    "\n",
    "## Objective and Prerequisites\n",
    "\n",
    "How can a mining company use mathematical optimization to identify which excavation locations to choose in order to maximize the gross margins of extracting ore? Try this modeling example to find out!\n",
    "\n",
    "This model is example 14 from the fifth edition of Model Building in Mathematical Programming by H. Paul Williams on pages 269-270 and 324-325.\n",
    "\n",
    "This example is at the intermediate level where we assume that you know Python and the Gurobi Python API and that you have some knowledge of building mathematical optimization models.\n",
    "\n",
    "**Download the Repository** <br /> \n",
    "You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem Description\n",
    "\n",
    "A company has obtained permission to conduct opencast mining within a square plot 200 ft $\\times$ 200 ft. The angle of slip of the soil is such that it is not possible for the sides of the excavation to be steeper than 45 degrees. The company has obtained estimates for the value of the ore in various places at various depths. Bearing in mind the restrictions imposed by the angle of slip, the company decides to consider the problem as one of the extracting of rectangular blocks. Each block has horizontal dimensions 50 ft $\\times$ 50 ft and a vertical dimension of 25 ft. If the blocks are chosen to lie above one another then it is only possible to excavate blocks forming an upturned pyramid. \n",
    "The three dimensional representation below shows four levels of excavation. We have numbered, in black, each block at each level, and the number in red represents the block underneath the four blocks of the level. For example, block 17 of level 2 lies underneath the blocks 1,2,5,and 6 of level 1.\n",
    "![pyramid](extractionPyramid.PNG)\n",
    "\n",
    "The profit for the extraction of ore at each block has been estimated. The goal is to find an ore extraction plan that maximizes total profit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Formulation\n",
    "\n",
    "### Sets and Indices\n",
    "\n",
    "$b,b2 \\in \\text{Blocks}=\\{1,...,30 \\}$.\n",
    "\n",
    "### Parameters\n",
    "\n",
    "$\\text{profit}_{b} \\in \\mathbb{R}^+$: Profit from extracting ore from block $b$.\n",
    "\n",
    "$(b,b2) \\in Arcs = Blocks \\times Blocks$: This parameter represent the arcs in the series-parallel graph describing the rules of extraction. The arc $(b,b2)$ in the adjacency matrix of this series-parallel graph has a value of 1 if block b2 is one of the four blocks above block b, and 0 otherwise. For example, arc $(29,24)$ represents that block 24 is one of the four blocks above block 29.\n",
    "\n",
    "### Decision Variables\n",
    "\n",
    "$\\text{extract}_{b} \\in \\{0,1\\}$: This binary variable is equal 1, if block $b$ is selected, and 0 otherwise.\n",
    "\n",
    "### Constraints\n",
    "\n",
    "**Extraction**: If a block is extracted, then the four blocks above it must also be extracted..\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{extract}_{b2} \\geq \\text{extract}_{b} \\quad \\forall (b,b2) \\in \\text{Arcs}\n",
    "\\end{equation}\n",
    "\n",
    "### Objective Function\n",
    "\n",
    "**Profits**: Maximize profits from the extraction of ore.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Maximize} \\quad \\sum_{b \\in Blocks} \\text{profit}_{b}*\\text{extract}_{b}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python Implementation\n",
    "\n",
    "We import the Gurobi Python Module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "\n",
    "# tested with Python 3.11 & Gurobi 11.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input data\n",
    "\n",
    "We define all the input data for the model and other Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a dictionary to capture the profit of the extracting of ore at each block.\n",
    "\n",
    "blocks, profit = gp.multidict({\n",
    "    ('1'): 0,\n",
    "    ('2'): 0,\n",
    "    ('3'): 0,\n",
    "    ('4'): -1500,\n",
    "    ('5'): 0,\n",
    "    ('6'): 1000,\n",
    "    ('7'): 0,\n",
    "    ('8'): -1500,\n",
    "    ('9'): -1000,\n",
    "    ('10'): -1000,\n",
    "    ('11'): -1500,\n",
    "    ('12'): -2000,\n",
    "    ('13'): -1500,\n",
    "    ('14'): -1500,\n",
    "    ('15'): -2000,\n",
    "    ('16'): -2500,\n",
    "    ('17'): 2000,\n",
    "    ('18'): 2000,\n",
    "    ('19'): -2000,\n",
    "    ('20'): 0,\n",
    "    ('21'): 0,\n",
    "    ('22'): -4000,\n",
    "    ('23'): -2000,\n",
    "    ('24'): -2000,\n",
    "    ('25'): -5000,\n",
    "    ('26'): 16000,\n",
    "    ('27'): 4000,\n",
    "    ('28'): 2000,\n",
    "    ('29'): 0,\n",
    "    ('30'): 2000\n",
    "})\n",
    "\n",
    "# Create a dictionary for the adjacency matrix of the series-parallel graph.\n",
    "\n",
    "arcs, value = gp.multidict({\n",
    "    ('30','26'): 1,\n",
    "    ('30','27'): 1,\n",
    "    ('30','28'): 1,\n",
    "    ('30','29'): 1,\n",
    "    ('29','21'): 1,\n",
    "    ('29','22'): 1,\n",
    "    ('29','24'): 1,\n",
    "    ('29','25'): 1,\n",
    "    ('28','20'): 1,\n",
    "    ('28','21'): 1,\n",
    "    ('28','23'): 1,\n",
    "    ('28','24'): 1,\n",
    "    ('27','18'): 1,\n",
    "    ('27','19'): 1,\n",
    "    ('27','21'): 1,\n",
    "    ('27','22'): 1,\n",
    "    ('26','17'): 1,\n",
    "    ('26','18'): 1,\n",
    "    ('26','20'): 1,\n",
    "    ('26','21'): 1,\n",
    "    ('25','11'): 1,\n",
    "    ('25','12'): 1,\n",
    "    ('25','15'): 1,\n",
    "    ('25','16'): 1,\n",
    "    ('24','10'): 1,\n",
    "    ('24','11'): 1,\n",
    "    ('24','14'): 1,\n",
    "    ('24','15'): 1,\n",
    "    ('23','9'): 1,\n",
    "    ('23','10'): 1,\n",
    "    ('23','13'): 1,\n",
    "    ('23','14'): 1,\n",
    "    ('22','7'): 1,\n",
    "    ('22','8'): 1,\n",
    "    ('22','11'): 1,\n",
    "    ('22','12'): 1,\n",
    "    ('21','6'): 1,\n",
    "    ('21','7'): 1,\n",
    "    ('21','10'): 1,\n",
    "    ('21','11'): 1,\n",
    "    ('20','5'): 1,\n",
    "    ('20','6'): 1,\n",
    "    ('20','9'): 1,\n",
    "    ('20','10'): 1,\n",
    "    ('19','3'): 1,\n",
    "    ('19','4'): 1,\n",
    "    ('19','7'): 1,\n",
    "    ('19','8'): 1,\n",
    "    ('18','2'): 1,\n",
    "    ('18','3'): 1,\n",
    "    ('18','6'): 1,\n",
    "    ('18','7'): 1,\n",
    "    ('17','1'): 1,\n",
    "    ('17','2'): 1,\n",
    "    ('17','5'): 1,\n",
    "    ('17','6'): 1\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Deployment\n",
    "\n",
    "We create a model and the variables. These binary decision variables define which block to extract ore from.\n",
    "\n",
    "Notice that the matrix of coefficients of the constraints is totally unimodular, therefore the decision variables can be defined in the interval $[0,1]$ and the problem can be solved as a linear programming problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file c:\\gurobi\\gurobi.lic\n"
     ]
    }
   ],
   "source": [
    "model = gp.Model('opencastMining')\n",
    "\n",
    "# Decision variable to extract ore from block\n",
    "extract = model.addVars(blocks, ub=1, vtype=GRB.CONTINUOUS, name=\"extract\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following constraints ensure that if a block is extracted, then the four blocks above it must also be extracted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extraction constraints\n",
    "\n",
    "extractionConstrs = model.addConstrs( (extract[b] <= extract[b2]  for b,b2 in arcs), name='extractionConstrs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to maximize the profits from the extraction of ore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Objective function\n",
    "\n",
    "extractionProfit = gp.quicksum(profit[b]*extract[b] for b in blocks )\n",
    "\n",
    "model.setObjective(extractionProfit, GRB.MAXIMIZE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 56 rows, 30 columns and 112 nonzeros\n",
      "Model fingerprint: 0x83427493\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [1e+03, 2e+04]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [0e+00, 0e+00]\n",
      "Presolve removed 22 rows and 12 columns\n",
      "Presolve time: 0.01s\n",
      "Presolved: 34 rows, 18 columns, 68 nonzeros\n",
      "\n",
      "Iteration    Objective       Primal Inf.    Dual Inf.      Time\n",
      "       0    2.3014500e+04   1.100300e+01   0.000000e+00      0s\n",
      "       9    1.7500000e+04   0.000000e+00   0.000000e+00      0s\n",
      "\n",
      "Solved in 9 iterations and 0.01 seconds\n",
      "Optimal objective  1.750000000e+04\n"
     ]
    }
   ],
   "source": [
    "# Verify model formulation\n",
    "\n",
    "model.write('opencastMining.lp')\n",
    "\n",
    "# Run optimization engine\n",
    "\n",
    "model.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "The total profit generated from the optimal ore extraction plan is $\\$17,500.00$. \n",
    "The block to extract ore from and its associated profit or loss are shown in the following table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Block</th>\n",
       "      <th>Profit/Loss</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>1</td>\n",
       "      <td>$0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>2</td>\n",
       "      <td>$0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>3</td>\n",
       "      <td>$0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>5</td>\n",
       "      <td>$0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>6</td>\n",
       "      <td>$1,000.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>7</td>\n",
       "      <td>$0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>9</td>\n",
       "      <td>$-1,000.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>10</td>\n",
       "      <td>$-1,000.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>11</td>\n",
       "      <td>$-1,500.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>17</td>\n",
       "      <td>$2,000.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>18</td>\n",
       "      <td>$2,000.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>20</td>\n",
       "      <td>$0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>21</td>\n",
       "      <td>$0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <td>26</td>\n",
       "      <td>$16,000.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       " Block Profit/Loss\n",
       "     1       $0.00\n",
       "     2       $0.00\n",
       "     3       $0.00\n",
       "     5       $0.00\n",
       "     6   $1,000.00\n",
       "     7       $0.00\n",
       "     9  $-1,000.00\n",
       "    10  $-1,000.00\n",
       "    11  $-1,500.00\n",
       "    17   $2,000.00\n",
       "    18   $2,000.00\n",
       "    20       $0.00\n",
       "    21       $0.00\n",
       "    26  $16,000.00"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Output reports\n",
    "\n",
    "extraction_plan = pd.DataFrame(\n",
    "    [(b, f\"${profit[b]*round(extract[b].x):,.2f}\") for b in blocks if (extract[b].x > 0.5)],\n",
    "    columns = [\"Block\", \"Profit/Loss\"],\n",
    ")\n",
    "extraction_plan.index=[''] * len(extraction_plan)\n",
    "extraction_plan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "H. Paul Williams, Model Building in Mathematical Programming, fifth edition.\n",
    "\n",
    "Copyright © 2020 Gurobi Optimization, LLC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
